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Abstract 


Flow and turbulence models applied to the problem of shock buffet 
onset are studied. The accuracy of the interactive boundary layer 
and the thin-layer Navier-Stokes equations solved with recent upwind 
techniques using similar transport field equation turbulence models is 
assessed for standard steady test cases , including conditions having 
significant shock separation. The two methods are found to compare 
well in the shock buffet onset region of a supercritical airfoil that 
involves strong trailing-edge separation. A computational analysis 
using the interactive boundary layer has revealed a Reynolds seal - 
ing effect in the shock buffet onset of the supercritical airfoil, which 
compares well with experiment . The methods are next applied to a 
conventional airfoil. Steady shock- separated computations of the con- 
ventional airfoil with the two methods compare well with experiment. 

Although the interactive boundary layer computations in the shock 
buffet region compare well with experiment for the conventional air- 
foil, the thin-layer Navier-Stokes computations do not. These find- 
ings are discussed in connection with possible mechanisms important 
in the onset of shock buffet and the constraints imposed by current 
numerical modeling techniques. 

Introduction 

Shock buffet or shock-induced oscillation (SIO) is large-scale flow-induced shock motion 
that involves alternating separation and reattachment of a boundary layer. In several recent 
computational studies, prominent features of the shock buffet of the 18-percent-thick circular- 
arc airfoil have been computed with Navier-Stokes and thin-layer Navier-Stokes codes (refs. 1 
and 2). Those studies highlighted the sensitivity of this problem to the type of turbulence 
and flow model and the importance of shock and trailing-edge separation in the onset of shock 
buffet. Although details of the shock buffet are sensitive to these factors, all computations 
have computed the onset Mach number for the circular-arc airfoil quite accurately. After the 
comprehensive time accurate calculations made for the shock buffet of the 18-percent circular-arc 
airfoil in reference 2, an assessment using current methods and turbulence models of predictive 
capabilities for several more widely used airfoils was undertaken. The present report shows the 
results of a computational study of this problem with both the interactive boundary layer (IBL) 
method and a thin-layer Navier-Stokes (TLNS) code. 

The physical mechanisms important in this problem can be investigated from a variety of 
viewpoints. For instance, shock strength is implicated in the identification of a Mach number 
range ahead of the shock for the 14-percent circular-arc airfoil in which shock buffet occurs 
(ref. 3). Geometry and trailing-edge viscous-inviscid interaction play a role as well. The 
18-percent circular-arc airfoil has trailing-edge separation prior to shock separation and shock 
buffet onset (ref. 4). Trailing-edge separation has long been associated with the onset of shock 
buffet. (See refs. 5 and 6.) Shock buffet for this airfoil is antisymmetric and displays hysteresis 
in the onset Mach number range, the latter of which is discussed in reference 7 in connection 
with the coalescing of a shock and trailing-edge separation. Questions remain, however, as to 
the important mechanisms involved for other airfoils. For instance, the NACA 0012 airfoil has a 
much weaker trailing-edge separation in the onset region and experiences one-sided shock buffet 
(ref. 8). Nor does onset for this airfoil have a hysteresis in Mach number, and it does not 
apparently display the sensitivity of the shock buffet range on Reynolds number that is evident 
for the 18-percent circular-arc airfoil (ref. 8). 



Experimental measurement of shock buffet onset is of course complicated by external effects 
such as wind tunnel noise, Reynolds scaling, and walls. However, experiment and several 
computations show transonic Mach numbers within an angle-of-attack envelope for the NACA 
0012 airfoil where shock motion inteasity and chordwise extent change from a localized shock 
oscillatory (or steady in the case of the computations) to a large-scale motion displaying limit 
cycle behavior. This has been studied experimentally for the NACA 0012 airfoil in reference 8, 
which represents an effort to provide quality steady and unsteady lifting surface results with 
minimal interference effects. In the test of reference 8, tunnel walls were contoured to match free 
air streamlines for nominal test conditions derived from Navier-Stokes computations. A much 
less extensive study of the NACA 00 12 airfoil has been presented in reference 9, which also reveals 
shock buffet behavior. Reference 10, in contrast, presents experimental data from a slotted- wall 
wind tunnel for the same airfoil and range of conditions that are steady. Experimental studies 
through the onset Mach number range for several supercritical airfoils confirm that these airfoils 
can also experience shock buffet (refe. 9, 11, 12, and 13). In summary, although difficulties 
remain in verifying onset and sorting out the various extraneous effects, it is clear that under 
the right conditions some conventional and supercritical airfoils experience shock buffet. 

The computations vary somewhat for airfoils other than the 18-percent, circular arc. Steady 
interactive boundary layer and Navier-Stokes solutions of the NACA 0012 airfoil have been 
previously published (e.g., refs. 14 and 15) and compared with the steady data of reference 10. 
Shock buffet interactive boundary layer computations for the same airfoil have been shown in 
previous publications with a time accurate integral boundary layer and the classical transonic 
small disturbance (TSD) equation (ref. 16) and a TSD using an Euler-like streamwise flux and 
a steady integral boundary layer (ref. 7). This last reference has identified the onset behavior 
for the NACA 0012 airfoil as a Hopf bifurcation point where the solution changes from an 
equilibrium point to a limit cycle solution. The critical point or onset location is the point 
having a zero amplitude limit cycle solution. Whether a supercritical airfoil behaves like a 
conventional airfoil or like the 18-percent, circular-arc airfoil in this and other respects, however, 
remains to be ascertained. But the fact that the interactive boundary layer approach has given 
shock buffet, onset for the NACA 0012 airfoil that compares well with the onset of reference 8 
(e.g., refs. 7 and 17) encourages one to pursue further investigation with this method. Although 
it is generally accepted that the boundary layer assumption is violated in many problems of 
this type, the interactive boundary layer method does make possible a broader study of the 
problem due to its efficiency. That is done here with a recently developed interactive boundary 
layer method using the CAP-TSD (Computational Aeroelasticity Program) potential code with 
a modified streamwise flux (ref. 18) and an unsteady compressible boundary layer solved in 
finite difference form. This method is shown to give very accurate results for many widely 
used attached and shock -separated steady test, conditions. Comparisons of wall shear, boundary 
layer velocity profiles, and pressure distributions are shown to match well with experiment and 
Navier-Stokes results. In view of the sensitivity of the 18-percent, circular-arc airfoil shock 
buffet to turbulence and flow model, comparisons of shock buffet onset for the NACA 0012 
airfoil using several turbulence models are shown. Results are presented for several variations 
of the k-u) turbulence model. The k-uj turbulence model embodies more flow 7 physics than one- 
or zero-equation turbulence models and is applicable to boundary layer dominated flows. It 
allows solution of the turbulence equations to the wall including the viscous sublayer and also 
allows modeling of free-stream turbulence and the effect of varying surface roughnesses. This 
allows the effect of these modeling parameters on shock buffet onset, to be investigated. The 
shear stress transport form of the model is used to compute details of the shock buffet of the 
NASA SC(2)-0714 airfoil (ref. 13). Comparisons with the experimental shock buffet data at. high 
Reynolds numbers of reference 13 are shown at several Reynolds numbers; this represents the 
first numerical study of the effect of turbulent boundary layer Reynolds number scaling on shock 
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buffet onset. The wind tunnel walls are not modeled computationally, and wind tunnel effects 
are only considered when using standard corrections to Mach number and angle of attack. 

Previous unsteady thin-layer Navier-Stokes shock buffet results for a conventional and 
supercritical airfoil using the Baldwin-Lomax model and flux differences using artificial smooth- 
ing were shown by other authors in references 9 and 19. In the present effort, the Navier-Stokes 
simulation of shock buffet flows about these airfoils is made with Roe’s upwind split flux differ- 
encing and an advanced turbulence model. Thin-layer Navier-Stokes computations are shown 
using the Menter k-u> shear stress transport (SST) and the Spalart-Allmaras (SA) turbulence 
models. Thin- layer Navier-Stokes results with a steady separated shock are shown that com- 
pare well with the steady interactive boundary layer results; these are followed with thin-layer 
Navier-Stokes computations in the shock buffet onset regions of the two airfoils. The effect of 
transition location in the shock buffet onset region of the NACA 0012 airfoil is also examined. 

First, the equations and some pertinent theory behind the methods and turbulence models 
used in this study are given. Then, the numerical method of the interactive boundary layer model 
is described. Last, computational results and comparisons with experiment are presented. 

Leading Order Matching and Boundary Layer Equations 

As the starting point for obtaining the boundary layer equations, the Navier-Stokes equations 
are nondimensionalized by the airfoil chord c* and the free-stream density, velocity, and viscosity, 
denoted by and , respectively. Nondimensional variables are defined by 
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= y l 

c* 



u* 

_ u* 


_ T* 


u — 


v ~ u*~ 

u oo 

T 

— y* 

x oo 



£_ 

p*-p% 



* 

p 

p — 

Poo 



P — 

Poo 


where u is nondimensional x- velocity component and v is nondimensional y- velocity component. 
Because of the dependence of the resulting equations on Reynolds number, Re = p^U^c* / 
the unsteady compressible boundary layer equations may be obtained upon expansion of 
quantities in terms of successively smaller orders of the parameter e. The boundary layer 
vertical scale is defined by y = cY , for a constant y as e — ► 0 and streamwise extent of 0(1). In 
the boundary layer, expansions for velocities, pressure, density, temperature, and displacement 
thickness are 

u ~ U\ + dJ 2 + - • ‘ v ~ eVi + + • • * 

p ~ Pi + eP-2 + * * • p ~ R\ + + * • • 

T ~ ©i 4* e©2 + * • • 6 ~ + • • • 

Displacement thickness is defined here as 



where the subscript e refers to values at the boundary layer edge. To match with the boundary 
layer, outer flow quantities are expanded as 


3 



u ~ u\ + CU '2 + * • • v ~ v\ + cv 2 + • • • 

p ~ Pi + cp 2 + * * * P ~ Pi + e P'2 + * • • 

T ~ Ti + eT'2 + ■ • • 

For the present purpose of obtaining equations usable in solving an airfoil problem, the boundary 
layer equations are Favre mass averaged and then Prandtl transposed. The coordinate and 
velocity definitions are Y = ±Y± + F±(x,t) and V\ = ±V + U\F± X + F±t for the upper and 
lower surface boundary layers, respectively. Here the airfoil surface height is defined as f± = cF±. 
The boundary layer equations remain unchanged in form, whereas the inviscid velocity injection 
after transformation has the changes to be noted subsequently. 

The equations are transformed using Levy-Lees variables. The transformation is 

f = I 

and ^ 

F = Ui V = -V^i§^(v^i m J Q ±Fdr ! 
p = P\ 6 = 0i 

where a/ = 2£. The boundary layer equations become 


V rf + F + Qf/Fjc — 0 

[(' + Wl , - - VF, - a,F, + FV, + F 2 = 0 

Pn =0 


( 1 ) 


In these equations l = pp and /j> = ppj*, where p is molecular viscosity and pj* is eddy viscosity. 


The complete set of adiabatic wall boundary conditions at p± =■ 0 are F = V = 0, p^± = 0, 
and 6 r) = 0. In the wake, the boundary condition is F r) ± = 0 for a symmetric wake at p± = 0, 
whereas an additional equation, such as ^-momentum, would be required for an asymmetric 
wake. 


The leading order matching conditions can be rewritten as 

F(Z>V±’i) = «i(*,0 ±,t) = 1+ <j> x (x,Q ±,t) (2) 

V - T) ± V r , = V ^7 ^ [peMh + $)] ( 3 ) 


as 7]± — ► oo and 

p(Z,t) =pi(x,0±,t) | 
0(£i»7±,<) = T\{x,0 ±,t) J 


(4) 


where ( u , v) = V<p. In equation (3), 9 = ( 1 — 0) dp. The leading order inviscid injection 

velocity is now 


tyy{x , 0^, t) 


j dipeUeS 0 
\ dx 



- P xt Or, 0 ,t)]dY 


Pe 1 + f±x + f±t 


(5) 
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as Y± — * oo. The plus and minus in equation (5) refer to the upper and lower surfaces, 
respectively. In the computations that follow we set f±i = 0. 

Boundary Layer Turbulence Model Equations 

Variables are nondimensionalized in a manner identical to the boundary layer equations, 
with turbulence kinetic energy and specific dissipation rate nondimensionalized by and 
U^/c*, respectively. The resulting equations are Favre mass averaged and the nondime nsional 
k-u> equations rescaled in the same way as the boundary layer equations with the largest order 
terms retained. The equations are transformed with a steady Prandtl transposition given in the 
previous section. 

Transformation with boundary layer variables defined in previous section brings the equations 
into a nonasymptotic nonsimilarity form that is compatible with the boundary layer equations. 
The form of the equations is 


dk_ 

dt 


h - ^ 


' . v dk ' 

( l+<T k,nb)frj 


du> o __ \dkdu d f du 

— - Pu. - a,/3nu + 2Re«7u,,„-— — F\ + — [(/ + <r w , nl T )-. ^ 


(6a) 


(6b) 


where the transformed convective derivative is defined d()/dt = a^d/dt + Fd/di)-k V d/dr}. 
The first two terms on the right, Pf. and P M , are transformed turbulence production due to 
mean flow gradients with Reynolds stresses replaced by using the Boussinesq approximation. 
The tensor form of the production terms (in terms of physical quantities) is given by 


Pk -**9^ 

2 

Rij — vtSij — Trkbij 

eu, 

Pu> = In-r 1 ; 

OX j 


( 7 ) 


where Sjj is the mean velocity rate of strain tensor and Rij is the modeled Reynolds stress. 
In this study, only the leading order boundary layer production terms have been retained. 
Furthermore, transformation of these equations results in velocity as well as density gradient 
terms arising from the normal Reynolds stress production term. The part of this production 
term involving density gradient has also been neglected to simplify the expression; this is justified 
because the effect of density gradients is expected to be relatively small in the present transonic 
applications. Also, in many models, the effect of the normal Reynolds stress production is 
neglected altogether, although in the present implementation the contribution due to velocity- 
gradient has been retained. This is based on the observation of Delery from experimental data: 
this contribution due to velocity gradient can become important in shock separating flows. (See 
ref. 20.) The remaining terms on the right of equations (6) result from diffusion and dissipation. 
After solution of the complete equation set, the eddy viscosity is found by vj = Re k/u>. 
Note that this equation set treats eddy viscosity as a scalar quantity and, therefore, does not 
correctly incorporate the relationship between ((t/) 2 ) and this can be important at a 

shock separation. 

Variations of the k-u> turbulence model can be created from equations (6) by using differ- 
ent sets of coefficients. The three variations of the k-ui turbulence model and the coefficients 
used in each is summarized in table 1. The first variation is the Wilcox k- ui model (with F\ in 
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Table 1. Coefficients for Boundary Layer Two-Equation Turbulence Model 


Model 

11 

a k\v 

( Tvj.il 

d„ 

k^ 

1 

0.5 

0.500 

0.750 

k-u )j k-t 

1 

0.5 

0,500 

0.750 


2 

1.0 

0.856 

0.828 

SST 

1 

0.85 

0.650 

0.750 


2 

1.00 

0.856 

0.828 


equations (6) set to 0). The second set is the k-u )/k-e model. In this variant the k-u model is 
used near the wall, whereas the k-e model is used in the outer boundary layer. The k-e model 
is activated by the factor F\ in equation (6b) (with F\ — 1). The coefficients of this two-layer 
variation are 


a kx = 0.5 

*k, 2= 

0 * = 0.09 


<7^1 = 0.500 
a u >,2 = 0-856 
K = 0.41 


01 = 0.0750 ' 

02 = 0.0828 


(8) 


7 n — 0n/0* ~ &u ! , n K? / \J 0* 


for n— 1,2. The k-u model (first variation) uses the first set of coefficients (n = 1). In the 
two-layer model ( k-u f k-e ), the coefficients of equations (8) are blended from layer 1 to 2 by 
using the function 


• • • (<Pn = 

The blending function is determined from F\ — \[l + tanh(argi)]. Menter and Rumsey give a 
function for arg^ involving variables y , k, and other constants. (See ref. 21.) In the present 
method a simple expression is used based on boundary layer thickness that maintains a nearly 
constant relative location (at approximately the center of the boundary layer) and width of the 
transitional region. The blending function for several representative velocity profiles is shown in 
reference 17. 

The third turbulence model is the SST model of Menter. (See ref. 22.) This model is based 
on the observation that in an adverse pressure gradient, the ratio of turbulence shear stress to 
kinetic energy is nearly constant. Although not applicable to wake reattachment, wall jets, or 
free turbulent flows, it should give good results in the present application. This model is again 
a two-layer model composed of equations (6) but with the coefficients for layers 1 and 2 now 
given by 

cr *. j = 0.85 i = 0.650 /?i = 0.0750 

<7*2 = 1-00 <7^2 = 0.856 02 = 0.0828 

The key difference of this SST model with the previous two models is its behavior in an adverse 
pressure gradient, where the model switches to a shear stress transport model with r = Gqfc; 
this is accomplished by defining eddy viscosity by 

aik Re 

i/ r f = 

max(aiu;, J F2^) 

with the constant a\ — 0.31 identical to the value used by Menter and Rumsey (ref. 21). The 
vorticity Q is given by du/dy. In the present formulation F 2 equals 1 throughout the boundary 
layer. 
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Wilcox in reference 23 shows through an asymptotic analysis of the dissipation rate equation 
that a balance of the diffusion and dissipation terms in the logarithmic region of an attached 
turbulent boundary as — ► 0 yields 


_ 61/ 

~ h (»+) 2 


This boundary condition is valid for the law of the wall and consistent with viscous sublayer 
behavior. In a separated turbulent boundary layer, the law of the wall is no longer uniformly 
correct, although the prevailing condition of a back flow region is one of turbulence interacting 
independent of mean flow (refs. 24 and 25). The same balance of terms gives the same result with 
k « (y + ) m for m required to satisfy outer back flow boundary conditions. The same expression 
is, thus, equally valid for attached and separated turbulent boundary layers. It should also model 
viscous sublayer behavior without resorting to implementation of a wall function or damping 
terms. In the present application, a variation of an approximation suggested in reference 26 is 
used, in which the dissipation rate at the airfoil surface is given by 

, 6O/1 . v 

(9) 

for normal mesh spacing A rj at the base of the boundary layer. Setting = 1 simulates the 
effect of a smooth surface, whereas smaller values simulate increasing roughness. 

At the outer edge of the boundary layer (77 — ► 00), a condition that, at leading order, k and u> 
remain finite requires that either their first and second derivatives in 77 tend to zero or that the 
two parameters be specified. One proposal (ref. 23) is to impose a zero gradient condition at 
the outer edge by solving the following equations: 


du) o 

F,^r=-l3’ku 


In this study the values of k and u are specified upstream and at the outer boimdary layer 
edge. This location makes nonzero turbulence energy and length scale gradients at the 
outer edge possible but does result in rapid convergence. One solution of these equations, 

uj = F e [/?£ 0 (£/£ 0 — l)]” 1 and k = C\F ^ 0 ^ ^(£/£ 0 — l) -/? with F e — Constant, can be used 
as an approximation for a more general case with £ 0 as an arbitrary constant of integration from 
which the initial value u?(£ a ) can be set. This solution represents a modeling of the decay of 
grid-generated turbulence of a flow moving at a uniform speed and therefore is also indirectly 
dependent on the upstream boundary of the solution domain. In an interacting boundary layer 
the problem is simplified. We can arbitrarily set £ 0 — 0 (i.e., airfoil leading edge ahead of the 
boundary layer) to ensure consistency of the results to follow for all grids. 

In a perturbation analysis of the defect layer equations, Wilcox (ref. 23) has k oc v^ T 
and lu oc u^./u e 6 in order to return an eddy viscosity proportional to u e 6. The dissipation and 
turbulence kinetic energy then are found by 


— j l r w/ p\ 


( 10 ) 


and k is derived by assuming an algebraic model value of eddy viscosity at the outer boundary 
layer edge. In equation (10), r w and u e 8 are implemented in a somewhat arbitrary way: 
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u t 8 = 5a/ 1 / 2 R e 1//2 and r w — 5a/ 1 / 2 Re 1 / 2 . The outer boundary value of kinematic eddy 
viscosity as used here is given by irp = (10) m z>TO ( m — “1> 0, 1) and = 0.089we^Re ly/2 . The 
parameter m allows varying the tree-stream kinematic eddy viscosity and kinetic energy. 


Thin-Layer Navier-Stokes Turbulence Model Equations 

Two turbulence models are used in the thin-layer Navier-Stokes computations. They are the 
Menter SST version of the k-uj turbulence model and the SA model. The manner in which the 
k-w model is implemented in the thin-layer Navier-Stokes code differs in several ways from the 
way it is implemented in the interactive boundary layer code. For that reason the k-u; model 
equations are repeated here. The Ar-u SST model equations can be written in the form 


Ok 


dk 


pw +pu jE~ = p k Re 1 -^^ Re + £r 


d 


dt 

du> 

dt 


dx 


dx 


J l 


die ' 

(M + *k,nlV)Tr- 

dxj A 


Re 


-1 


duo duu i 9 1 dk du) i 

p— + puj ^7- = Put Re 1 - f3puj z Re + 2 ^ — F\ Re 1 


dxj dxj 


+ 


dx: 


did 

(v + ct^j 


Re 


-1 


(ii) 


The eddy viscosity is given by 


. ( pk a \pk \ 

' ,3 ' =m,n b'w? R *J 

The production terms are approximated by 

Pk = PT& 1 
Pul = 1P& 1 

where Q is the magnitude of boundary layer vorticity du/dy. Note that this differs from 
the production terms used in the interactive boundary layer method where an additional 
production/destruction term is included. The effect of the difference would be to reduce the 
turbulence production in a favorable pressure gradient in the interactive boundary layer solution 
relative to that of the thin-layer Navier-Stokes solution. The constants are calculated from 
<j> = (1 — + Fi<j> 2, where <j>i and are the constants given here: 


°k, 1 = 

0.85 

<7^ 1 = 0.500 
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0.0750 

> 
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0.09 

F\ = 1 — tanh(T 4 ) 
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Note that one of the coefficients and the blending function used in this implementation differ 
from that used in the interactive boundary layer implementation. 

The SA turbulence model equation is 


dd dv 

~di + u j~ = Cbl ^ ~ 


dxj 
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The following definitions are used in tins equation: 


(13) 


where 


/<2 = C,z exp( -C/4X 2 ) f w = 9 
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Eddy viscosity is calculated from the equation 

IH = pvf v l 


(14) 


(15) 


where 


f v 1 — 


_ v 
^ v 


X 3 +Ol 3 

The constants are 

C bl = 0. 1355 a = 2/3 C 6 2 = 0.622 k = 0.41 C w2 = 0.3 ) 

C w 2=2.0 0,1 = 7.1 03 = 1.2 04=0.5 C v 2 = 5.0 

Oi 1 + O2 

c w 1 1 


(16) 


Additional details on the thin-layer Navier-Stokes equations and the numerical method used can 
be found in reference 2. 

Numerical Method for Interactive Boundary Layer Equations 

The quasi-simultaneous interactive boundary layer method as implemented by Davis and 
Werle has been modified for the present application. It couples the unsteady boundary layer 
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equations solved in finite difference form and the CAP-TSD potential code via a consistent 
matching condition and interaction law. The TSD solution uses a streamwise flux that 
approximates an Euler flux to several orders. The turbulence equations are decoupled from 
the boundary layer equations and solved first at each time step n + 1. Mean flow values at the 
previous time step n are used in the turbulence model computation. The equations are solved 
with full implicit finite differencing. Both the streamwise and normal derivatives have upwind 
differencing to enhance diagonal dominance, and the nonlinear coefficients involving k and u are 
lagged at the previous time step, giving the turbulence modeling first-order accuracy in time 
and space. 

Mean flow quantities used in the turbulence equations and the values of k and u actually 
used in computing eddy viscosity have also been found to result in a more stable boundary 
layer interaction if smoothing is applied. An averaging function of the form <I> = £<*<£ is used 
with summation over 6 to 9 adjacent grids. In computing equations (6), <£ = [F, V] is replaced, 
whereas in computing eddy viscosity, $ = [k,tj] is replaced with averaged values of Note 
that in neither of these cases is the artificial smoothing applied to either the boundary layer or 
turbulence field solution directly; thus, the form of the leading order equations (eqs. (6)) being 
solved is maintained. 

For the computation of the finite difference equations, the following terms are defined: 


/ 
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The mesh spacing functions are defined by 

Aii = it - it - 1 a m = Vi+i- m 


and 


Ar?ave = | (At y + A rjj-i) 

The boundary layer equations are finite differenced as follows: 
Continuity: 



F 


n+1 
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X-momentum: 


liMlViS - ( £ «+l Mi + Vl/A + s 3i - 2AF”+>)F5« + Vi/A^l 

+ (?/{^+i + JjVX* 1 + gj VJil) - V ; f 1 [W 2i (FJ+ 1 1 - }*+>) 
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+ AF" +1 F"+‘ + S: U F[; + (1 - A)3 2 (?it‘ - F," +1 ) 

Y-momentum: 

Py+1 ~ Pi j — 0 

where L — l + Irp. In the X-momentum equation (and the turbulence model equations to be 
shown next) 

f F^. +1 


and 


A = max 


A = max 


i^T’ 0 


yn+i 

—! Li — 0 

\Vij +l \ . 


The tilde represents the value from the previous iteration. The matching condition is differenced 
as follows: 


yn+1 ^max / un+1 yn+1 

v ij max \ v ij max v ij max — 1 
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The term 6 is explicitly integrated through the boundary layer and is treated as a known quantity 
at each inversion. The interaction law is of the form 

u" +1 = Ci + Di (/>e«e<5” +1 ) 

where C and D result from the in viscid governing equation. Streamwise sweeps of the bound- 
ary layer are performed. At the outer edge, the X-momentum equation is replaced with the 
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interaction law and the Y-momentum equation replaced with the X-momentum equation. The 
turbulence model equations are differenced as follows: 


+ (A^ +1/2 fy + + **3/ + a li.0 

+ - *?;' ) + < i - -'w,,(-:; +i - w? +_>,)] 

+ Fg[Ao 1; ^f 1 - uf+i) 4 ( 1 - A )S 2i « + i- - u't ' )] = 

and 

- + <%♦«.*> + S-./» 7 ' + 3si + “«*W - V./» 7 '*!ft 

4 - *,?■)+ (i - WnK* 1 - ‘5-i)i 

4 ^[AS 1( (tj +1 - *”:") 4(1- AJSaftf+l - *;'+>)] = h 4 5,,-tg 

where Lq = / + The parameters A and A are as defined earlier for the boundary layer 

equations. These equations are replaced at the outer boundary layer edge with the boundary 
conditions discussed in a previous section. Solved in this form, the equations decouple resulting 
in two efficient scalar tridiagonal inversions through the boundary layer. 

One or two global subiterations of the turbulence equations at each time step typically are 
required to reach single precision machine accuracy. The complete code has been vectorized 
with an average performance on a CRAY Y-MP computer of 41 Mflops, with the in vise id 
solver requiring approximately 150 Mflops and the boundary layer solver requiring approximately 
25 Mflops. Run times for converged steady state solutions are summarized in table 2. They 
compare quite well with other viscous codes. 


Table 2. CPI Execution Times of IBL Code 


Grid size 

Equivalent 

C-grid 

min 

At 

Computer time, min 
(a) 

Coarse 

134 by 120 

0.03 

0.015 

5 

Medium 

2G0 by 120 

0.015 

0.008 

20 

Fine 

360 by 125 

0.008 

0.004 

55 


a To converge to within 2 percent, of c/. 


Steady Computations 
RAE 2822 Airfoil 

The accuracy of the interactive boundary layer model is assessed by comparing steady state 
computations with the data of reference 27 and integral interactive boundary layer and Navier- 
Stokes results. The k-cu SST turbulence model is used in both the present interactive boundary 
layer and thin-layer Navier-Stokes computations. The experimental database for Computer 
Program Assessment (ref. 27) provides steady pressure, boundary layer, and wake data for sev- 
eral transonic cases. Those to be calculated here are AGARD cases 6, 9, and 10, which range 
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from the relatively easy case to the fairly difficult case 10 that has quite a strong shock boundary 
layer interaction. This test also has well-known problems related to tunnel interference, which 
make corrections to the data necessary. Consequently, a degree of uncertainty exists, especially 
in case 10. 

Experimental transition was fixed at 3 percent chord. In the present computations, transition 
has been simulated by starting the boundary layer just behind the leading edge with a very 
small upstream eddy viscosity (or k in the k-u) model) and placing computational wall surface 
roughness, through boundary conditions in the turbulence model, from 3 to 5 percent chord. 

The first calculated results are for AGARD case 9 shown in figure 1. Except at the upper 
surface leading edge where the superiority of the Euler solution on a body-fitted coordinate 
system over the present TSD solution is evident, the present C p distribution is appreciably 
better than that of reference 28 using an algebraic turbulence model. The shock strength and 
location using the present interactive boundary layer method match experiment very well, where 
algebraic models typically overpredict both shock strength and location. This trend is also 
evident in the skin friction and displacement thickness distributions presented in figures 2 and 3. 
The present results in these figures match experiment better except possibly near the leading 
edge. Otherwise, the recovery aft of the shock is still slightly too strong. 

With this case, we can also make a comparison of accuracy and efficiency for grids having 
differing amounts of grid stretching in the boundary layer. Figure 4 presents c p distributions for 
boundary layers having an average of approximately 175 and 85 boundary layer grids normal to 
the airfoil surface. As shown in table 2, the coarser grid reaches moderate engineering accuracy 
fairly efficiently. For the sake of accuracy, all the previous and remaining results are with the 
finer boundary layer grid spacing. 



Figure 1. Pressure coefficients for RAE 2822 airfoil, AGARD case 9, at A/exp = 0.73, a = 3.19°, and Re = 6.5 X 10 6 . 
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Figure 2. Friction coefficients for RAE 2822 airfoil, AGARD case 9, at M exp = 0.73, a ~ 3.19°, and Re = 6 



Re = 6.5 x 10 6 . 



.5 x 10 6 . 


19°, and 


Figure 4. Effect of boundary layer grid stretching for AGARD case 9 at M eX n = 0.73, a = 3.19°, and Re = 6.5 x 10 ( 


Results using the present interactive boundary layer method are shown with experiment and 
with other computational results for AGARD case 6. As seen in figure 5, the present results 
are not quite as good as the very good match of the previous case, although the displacement 
thickness, seen in figure 6, appears to match experiment well. Since the present method employs 
a finite differencing of the boundary layer, velocity profiles are also presented in figure 7. Again, 
except where boundary layer recovery is somewhat too rapid, the results are very good. 



Figure 5. Pressure coefficients for RAE 2822 airfoi, AGARD case 6, at M exp = 0.725, a = 2.92°, and Re = 6.5 x 10 6 . 



Figure 6. Displacement thicknesses for RAE 2822 airfoil, AGARD case 6, at Mexp = 0.725, a = 2.92°, and Re =s 6.5 x 10 6 . 
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O Experiment 





The conditions so far, with mild viscous-inviscid interaction, have been quite easy. AGARD 
case 10 is more challenging and takes 3-4 times longer to compute than the previous cases. The 
integral interactive boundary layer computations of references 15 and 29 are presented in figure 8, 
and the Navier-Stokes computations of reference 30 using the Johnson-King model are also 
presented in figure 8. The computations of the other references are at M = 0.75 and a = 2.81°. 
With the present method, at a — 2.81°, the solution exhibited growing shock oscillations and 
required a reduction to a = 2.70° to reach a steady state. At this angle, the present results 
represent a very good match with both experiment and the Navier-Stokes results of reference 30. 
Finally, note that the accuracy of the present interactive boundary layer computations for all 
these AGARD cases, which is superior to the integral boundary layer results shown here, may 
sufficiently justify the present use of a finite differenced boundary layer rather than the empirical 
and sometimes ad hoc closure relations used in integral boundary layer methods. 
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Figure 8. Pressure coefficients for RAE 2822 airfoil, AGARD case 10, at M exp = 0.75 and Re = 6.2 x 10®. 

NACA 64 AO 10 Airfoil 

A steady state solution using the present interactive boundary layer method and the k-u 
SST model has been computed for the NACA 64A010 airfoil at an angle of attack of 2° 
Moo = 0.80, and Re = 2 x 10 6 . At an angle of attack of 2°, the computed flow field separates 
over approximately 6 percent of the chord, representing a case for which the boundary layer 
has significant shock separation and yet the available experimental data appear to be steady. 
Thin layer Navier-Stokes results from using the k-u SST turbulence model and experiment, 
both from reference 21, are presented in figure 9 for comparison with the presently computed 
interactive boundary layer results. The present SST model results compare quite well with the 
Navier-Stokes results except near the leading edge and at the shock location. 
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Figure 9. Pressure coefficients for NACA 64A010 airfoil at M* = 0.80, t* = 2.0°, and Re = 2 x 10 6 . 


NACA 0012 Airfoil 

The present interactive boundary layer and thin-layer Navier-Stokes computations both, 
using the k-u> SST turbulence model, are presented for a steady condition of reference 8 near 
shock buffet onset. The thin-layer Navier-Stokes computations use Roe's flux difference splitting 
with up-wind biased third-order differencing and flux limiter. Multigrid is used to accelerate 
convergence at each time step. The Ar-u; SST turbulence model used in the thin-layer Navier- 
Stokes computations is identical to that of reference 21. The only structural differences between 
that and the present interactive boundary layer implementation is the inclusion ha the interactive 
boundary layer model of the production term proportional to kinetic energy, slight differences 
in the coefficients, and a different blending function. 

For the Navier-Stokes computations, several C-grids were tried having various spacings. Each 
grid gave qualitatively similar results. The grid used in computing the results shown here is 
moderately fine, with dimensions of 297 by 121. This grid extends 10 chords downstream and 
9 chords away from the airfoil. Outer boundary conditions are for free air. Normal wall spacing 
is A & 5. The inner portion of this grid is shown in figure 10. The grid used in the interactive 
boundary layer computation is the fine grid discussed earlier. 

For reference purposes, the thin-layer Navier-Stokes pressure distribution and that from the 
present interactive boundary layer method are presented in figure 11 for a steady condition just 
below onset at M lX > = 0.775, a = 2.05°, and Re = 10 x 10 () . From 5 to 10 percent chord on 
the upper surface, aft of the shock, and a portion of the lower surface are the only areas where 
the two computations differ slightly. Since the authors of reference 2 found that turbulence has 
a significant impact on buffet computations, an examination of Reynolds stress levels for this 
steady condition by the two methods is in order. A comparison of the data in figure 12 shows 
that the interactive boundary layer gives peak values that are less than half those produced 
by the thin-layer Navier-Stokes implementation up to at least 5 percent chord. The thin-layer 
Navier-Stokes method on the other hand is producing relatively constant peak turbulence levels 
from the leading edge with no evidence of transition. The difference in levels at x = 0.44 reflects 
a slight difference in shock location. Otherwise a significantly higher turbulence level Is produced 
by the interactive boundary layer in the outer wake area of the boundary layer combined with 
a correspondingly thicker boundary layer from the shock to the trailing edge. 
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This case, just below the experimental shock buffet onset of reference 8, is presented in 
anticipation of the unsteady computations shown next. Although the thin-layer Navier-Stokes 
and interactive boundary layer computations of the shock buffet for the supercritical airfoil show 
good agreement, the shock buffet computations for the NACA 0012 airfoil with the two methods 
are shown next to give considerably different results. 



Figure 10. Near field of thin- layer Navier-Stokes grid. 



Figure 11. Pressure coefficients for NACA 0012 airfoil at M x . = 0.775, a = 2.05°, and Re = 10 x 10 6 . 
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Unsteady Computations 

NASA SC(2)-0714 Airfoil 

The present interactive boundary layer and thin-layer Navier-Stokes methods have been used 
to compute in the shock buffet onset region of the NASA SC(2)-0714 airfoil, corresponding to 
M ex p = 0.74 of figure 13. The conditions shown in that figure are from the experimental data 
set discussed in reference 13 from a high Reynolds number wind tunnel test conducted in the 
Langley 0.3-Meter Transonic Cryogenic Tunnel. The purpose of this test was to provide data 
of transonic conditions of a fixed and pitching NASA SC(2)-0714 airfoil through a range of 
Reynolds numbers. Reference 13 documents cases for fixed angles of attack from that test, some 
of which show shock buffet. These data are significant in the present context not only because 
shock buffet is displayed but because a slight Reynolds scaling is revealed in the onset location. 
The data require, however, significant correction for wind tunnel effects such as downwash or 
Mach number. For instance, based on the theory of reference 31, at c/ = 0.93, the angle of 
attack corrected to account for induced downwash would be A a = —1.6°. The corrected Mach 
numbers for this wind tunnel are found in reference 32. Only the Mach number corrections are 
used in the computations to follow. 

The experimental conditions showing shock buffet are cases 5 and 7 in figure 13 which are at 
M exp = 0.74, q = 3.0°, and Re= 15 x 10 C and 30 x 10 () , respectively. Those experimental data 
display a Reynolds scaling effect in the location of the shock buffet onset, a feature that is also 
shown in the present computations. In the computations that follow 7 , the critical point at which 
onset occurred was found by increasing the angle of attack by increments of 0.2° until onset 
occurs and interpolating to the bifurcation point. The bifurcation point location was inferred 
by the relative decay or growth rates of the shock oscillation at successive angles of attack. To 
assess the sensitivity of these computations to the grid used, this procedure w^as repeated with 
the three grids discussed previously. The buffet boundary was identical for all to within ±0.1°. 
Based on this result, the shock buffet interactive boundary layer computations were made wdth 
the finest of the three interactive boundary layer grids. 


O Cases 1-3; Re = 6 x 10 6 , 15 x 10 6 , and 30 x 10 6 
D Cases 4 and 6; Re = 15 x 10^ and 30 x 10^ 

^ Cases 5 and 7; Re = 15 x 10^ and 30 x 10^ 
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Figure 13. Experimental data (ret. 13) lor NASA SC(2)-0714 airfoil, steady and unsteady cases. 
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Interactive boundary layer results shown in the figures discussed subsequently are compared 
with experimental data from reference 13. As seen in figure 14, computed onset is very near 
the angle of attack of the experimental shock buffet at a = 3.0°, but the solution reveals a 
slight Reynolds number effect. This effect is made clearer by the unsteady pressure coefficient 
moduli plotted in figure 15 for solutions at successive angles of attack. Quite near buffet onset, a 
difference in intensity at different Reynolds numbers exists whereas deeper into buffet, difference 
diminishes. The similarity of the solutions near buffet onset with the data shown in figure 16 
suggests that the experimental conditions are also very near shock buffet onset. Figure 17 
presents computed and experimental upper surface mean pressure coefficient distributions for 
cases 5 and 7. The computed mean c p chordwise distributions in figure 17 and the normalized 
unsteady \c p \ distributions in figure 15 suggest that the angle of attack to match the experiment 
at a exp = 3.0° should be computed a.t an angle of attack close to that angle. On the other 
hand, figure 18 shows the computed frequency approaching experiment slightly as the angle of 
attack is increased to 3.4°. This feature and the fact that the corrected angle of attack based on 
reference 31 would be a C orr = 1.4° is somewhat puzzling. The experimentally observed rearward 
shift of the shock oscillations with increasing Reynolds number is also not represented in the 
numerical solutions. The computations place the region of shock oscillation several percent 
forward at the higher Reynolds number. 

Thin-layer Navier-Stokes computations with the SA model are shown next, for conditions near 
the shock buffet case 5 of figure 13 at a Reynolds number of 15 x 10^. The grid used had mesh 
spacing essentially the same as that used in the steady NACA 0012 thin-layer Navier-Stokes 
computations discussed earlier. Wall spacing is A y+ % 6. To initiate the thin-layer Navier- 
Stokes computations, a steady solut ion at an angle of attack several degrees below experimental 
buffet onset was obtained. Second-order time accurate computations were begun with step jumps 


O Present IBL, steady 

□ Present IBL, steady 
cases 4 and 6 


^ Present IBL, computed 
buffet onset 



Reynolds number 


Figure 14. Computed shock buffet ousel for NASA SC( 2 )- 07 14 airfoil, k- w SST model, at M com p = 0.725 and 
4 /ox p — 0. / 4. 
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in the angle of attack with an amplitude of ±2° to approximate a single shock oscillation; this 
was used to initiate the shock buffet. A time step size of 0.04 (nondimensionalized by speed of 
sound) was used for the time accurate computations with 4-5 subiterations per time step and 
multigrid. The computations were continued to either a steady state or a converged limit cycle 
oscillation. 


a comp» 

deg 



Figure 15. Shock buffeting normalized moduli for NASA SC(2)-0714 airfoil at M CO mp = 0.725 on upper surface. 
Present IBL. 



Figure 16. Shock buffet moduli for NASA SC(2)-0714 airfoil at peak frequencies for cases 5 and 7 at M^ p = 0.74 
and c*exp = 3.0° on upper surface. Experimental data from reference 13. 
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Figure 17. Shock buffet mean pressure coefficients for NASA SC(2)-0714 airfoil at Afcomp = 0.725 on upper surface. 

This procedure was repeated at successive angles of attack at increments of 1/2°. The results 
of these computations are shown in figures 18 to 22. These results compare reasonably well 
with experiment as well as with the IBL computations. Onset of shock buffet occurs between 
the angles of attack from 3° to 3.5°; this is compared with an onset at an angle of attack of 
2.97° from the IBL computation. The c/ trace of figure 20 shows shock oscillations at a = 3.5° 
having reached converged limit cycle behavior, the frequency of which is quite close to that of 
experiment. (See fig. 18.) The interactive boundary* layer and thin-layer Navier-Stokes mean 
Cp distributions in figure 21 are in moderate agreement, although the thin-layer Navier-Stokes 
computation is in somewhat better agreement with experiment in the shock region. The pressure 
recovery toward the trailing edge is also somewhat more rapid in the interactive boundary layer 
solution. Figure 22 shows experimental chordwise distributions of |?p| normalized to peak values 
of the fundamental along with those of the interactive boundary layer and thin-layer Navier- 
Stokes computations. The relative magnitudes of the computed fundamental and first harmonic 
match experiment very well, although in both the chordwise extent of the shock oscillation is 
forward of experiment. Both underestimate the overall extent as well. For example, the shock 
motion extent in the thin-layer Navier-Stokes computation is roughly half that of the experiment. 

This comparison of both computed results with experiment offers hope that both the 
interactive boundary layer and the thin-layer Navier-Stokes methods give a reasonably accurate 
computation of shock buffet onset. As seen in the next section this view may be premature. 
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O Present IBL (k - co SST) 
□ Present TLNS (S A) 

<3> Experiment (ref. 13) 


Figure 18. Shock buffet frequencies for NASA SC(2)-0714 airfoil for case 5. 

O Present TLNS, steady case 5 
□ Present TLNS, steady case 4 
+ Present TLNS, buffet 
4 r 

♦ Estimated experimental 
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Reynolds number 

Figure 19. Computed shock buffet onset for NASA SC(*2)-0714 airfoil, SA model, at A/comp = 0.725 and M 


c l 1.00 



Figure 20. Shock buffet lift coefficients for NASA SC(2)-0714 airfoil, SA model, at A/ ronip = 0.725, a co 
and Re = 15 x 10 6 . Present TLNS. 




NACA 0012 Airfoil 

The present interactive boundary layer and thin-layer Navier-Stokes methods have been used 
to compute in the shock buffet onset region of the NACA 0012 airfoil. In view of the sensitivity of 
shock buffet to turbulence modeling observed in other studies and the extensive data available 
for this airfoil, an analysis has been made of the effect of turbulence model parameters on 
buffet onset with the interactive boundary layer model. For the interactive boundary layer 
computations, the medium resolution grid was used in comparisons of turbulence models to 
reduce computation times. The time step was 0.012 with 2 to 3 boundary layer subiterations 
per time step. In assessing onset location for each of the turbulence models, the angle of attack 
was successively increased by increments of 0.2°. Onset was taken as the first angle at which 
computed oscillations continued to grow after several initial cycles of transient shock oscillations 
had passed. 

The buffet boundaries found by using the present interactive boundary layer method and 
other published results along with the experimental values from reference 8 are shown in figure 23. 
The k-u SST turbulence model uniformly gives the boundary lower than the other turbulence 
models and overall represents the best match of all the models. Both the Wilcox one and two 
layer models appear to give only slight improvement over results with an algebraic model. From 
figure 24, there is also a rather modest decrease in onset angle with decrease in the value of 
turbulence kinetic energy k at the boundary layer edge corresponding to a lower free-stream 
turbulence level. 

The thin-layer Navier-Stokes shock buffet onset computations for this airfoil were done with 
both the k-u SST and the SA turbulence models. The same procedure was used for these com- 
putations as for the supercritical airfoil. To initialize, a steady state solution was obtained at 
an angle of attack of 2° below experimental buffet onset. A steady lift coefficient was obtained. 
Second-order time accurate computations were begun by using time steps of 0.005 and 0.04 
with the k-u SST and SA turbulence models, respectively. Step jumps in the angle of attack 
with an amplitude of ±2° simulating a single shock oscillation were used to start the shock buffet. 


a, 

deg 


O fc-co SST steady 

• k - co SST unsteady 

■ k - co 2-layer Wilcox unsteady 

♦ k - co 1 -layer Wilcox unsteady 
▲ Algebraic, unsteady (ref. 16) 
^ Green's lag, unsteady (ref. 7) 




Figure 23. Computed shock buffet onset found by 
various turbulence models and experiment for 
NACA 0012 airfoil. IBL. 


F igure 24. Effect of boundary layer edge turbulence 
kinetic energy on shock buffet onset. Present 
IBL. 
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Computations were completed with the two turbulence models at the highest Mach number at 
a final angle of attack well into the buffet region. Because of a time step restriction, apparently 
imposed by the k-ui SST model, this is the only unsteady computation that has been made by 
using that model. Although comparison of the S A and the k-Lj SST model results are not shown, 
computation confirmed that the two models give very similar results in the shock buffet region 
and that the ensuing thin-layer Navier-Stokes computations could be concluded with the SA 
model alone. Accordingly, all unsteady results shown have been computed with the SA model. 

The results of this investigation are shown in figures 25 and 26. At M 0 0 = 0.775 and 
q = 4.0°, the ci oscillation amplitude had damped to about 10 percent of the peak excursion 
1 1/2 cycles after the driving force was eliminated. (See fig. 26(a).) Computations several cycles 
farther showed oscillations continuing to damp out. Although skin friction is not shown, the 
final steady state solution has a significant shock separation that marginally reattaches just 
before the trailing edge. This procedure was repeated at higher angles of attack with a similar 
level of damping. Computations at M 0 c = 0.75 over a range of angle of attack give slightly less 
damped shock oscillations, whereas at a Mach number of 0.725 the shock oscillations are even 
less damped. These results have been obtained with the turbulence model computed throughout 
the flow field. The effect of a turbulent transition occurring in the 10-percent-chord range has 
been simulated as well by turning the turbulence production terms on at 10 percent chord with 
zero turbulence production ahead of that point. This effectively creates a laminar flow in the 
leading edge region. As seen in figure 26(b), computing the turbulence from the leading edge 
and starting at 10 percent chord appeared to have little effect on the outcome. A variation in 
the level of damping is also seen at the lower Mach number. At a = 6°, just beyond onset, 
the oscillations die out slowly, whereas at a = 9° the oscillations are again strongly damped. 
(Compare fig. 26(b) with fig. 26(c).) It is interesting to note a variation in amplitude in the 
limit cycle shock buffet with angle of attack in the interactive boundary layer computations of 
reference 7 for this airfoil. 

In summary, thin-layer Navier-Stokes computations with the NASA SC(2)-0714 airfoil have 
shown sustained shock buffet that compares well with experiment. The present thin-layer Navier- 
Stokes computations with the NACA 0012 airfoil have shown moderate to strongly damped shock 



Figure 25. Shock buffet onset for NACA 0012 airfoil. 
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oscillation, with no evidence of a developing limit cycle shock buffet. Turbulence transition 
location does not significantly alter this situation for the NACA 0012 airfoil. A limited effort 
at refinement of the TLNS grid likewise has not resulted in sustained shock oscillations. As 
for differences in the solutions of the two airfoils, clearly the supercritical airfoil has a stronger 
viscous-inviscid interaction behind the shock than does the conventional airfoil; this can be seen 
in the comparison of the skin friction for the two airfoils shown in figure 27. This and the fact 
that differences in the amount of trailing edge separation due to Reynolds scaling effects result 
in somewhat different onset locations at different Reynolds numbers for the NASA SC(2)-0714 
airfoil tend to confirm that trailing edge interaction has a considerable effect on the onset of 
shock buffet. Geometry certainly plays a role, the effect of which has not been fully explored in 
this report. 



Figure 27. Skin friction coefficients for two airfoils at M. x . = 0.725 and tv = 3.0°. Present TLNS. 


Concluding Remarks 

An interactive boundary layer method using a transonic small disturbance, potential outer 
flow model with an Euler-like streamwise flux has been coupled with several variations of the 
k-w turbulence model. This method has been found to compute very accurately many standard 
steady transonic test cases. Several steady computations have shown that the interactive 
boundary layer method is capable of giving results that compare very well with thin-layer Navier- 
Stokes results. Turbulence levels differed in some minor respects between the two methods, but 
in general, these results confirm that the interactive boundary layer method is capable of quite 
good accuracy. This accuracy, which is superior to many other integral boundary layer results, 
is also justification for the use of a finite differenced boundary layer, showing the boundary layer 
method at its best, rather than resorting to the more widely used empirical and sometimes ad 
hoc closure relations used in integral boundary layer methods. 

With the interactive boundary layer method, a study has been made of the flow and 
turbulence modeling necessary to accurately model shock buffet onset. Both the interactive 
boundary layer and a thin-layer Navier-Stokes method have been employed. Computations using 
the two methods in the shock buffet region of a 14- percent-thick supercritical airfoil compare 
well with experiment. These results suggest that both methods are capable of modeling shock 
buffet onset of that airfoil quite well. 
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Variants of the k-cu turbulence model with the interactive boundary layer method have been 
used to compute shock buffet onset for the NACA 0012 airfoil. These solutions confirm that 
turbulence model has an influence on the accuracy of the computed onset location most notably 
at higher Mach numbers. When comparing the interactive boundary layer and thin-layer Navier- 
Stokes computations, the results from the two methods were found to be quite different when 
computing the shock buffet of the NACA 0012 airfoil. The shock buffet onset computed with the 
interactive boundary layer and the k-u) SST (shear stress transport) turbulence model compares 
well with the onset of shock buffet seen in the data of NASA TP-2485. In contrast, the present 
thin-laver Navier-Stokes computations have uniformly shown damped shock oscillations well 
into the shock buffet region. Transition locations fixed respectively at the leading edge and at 
10 percent chord w 7 ere found to yield qualitatively similar thin-laver Navier-Stokes results. 

From the computed results, the supercritical airfoil has a much stronger viscous-inviscid 
interaction behind the shock than the conventional airfoil. Trailing- edge viscous-inviscid 
interaction has been shown by these results to have a considerable effect on the onset of shock 
buffet. What remains to be assessed computationally is the influence of wind tunnel walls and 
the effect of the numerical accuracy of the thin-layer Navier-Stokes scheme used in this report. 


NASA Langley Research Center 
Hampton, VA 23681-2199 
October 8, 1997 
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